      SUBROUTINE NEWMAG
     &          (spheroid,alt,rlng,rclt,bmf,dip,b,br,bp,bt)

c***********************************************************************
c                   subroutine newmag
c***********************************************************************

c  Program Source:  Naval Ocean Systems Center - Code 542

c  Date:
c     01 May 1991

c  Function:
c     Calculates the geomagnetic field values from a spherical harmonic
c     model. This routine was developed following the form found in the
c     routine named NEWMAG written by McIlwain of the University of
c     California at San Diego. This new program is designated NEWMAG85
c     but uses the same subroutine name and argument list as NEWMAG for
c     compatibility. It incorporates the spheroidal earth correction
c     used by the the World Data Center's collection of routines which
c     we have designated GEOMAG. These routines were written by

c        Bill Flanagan
c        National Geophysical Data Center
c        World Data Center-A for Solid Earth Geophysics
c        NOAA, E/GC1, 325 Broadway
c        Boulder, CO  80303

c     and

c        Norman W. Peddie, U.S. Geological Survey, Mail Stop 964,
c        Federal Center, Box 25046, Denver, CO 80225

c        who based his routine on 'IGRF' by D. R. Barraclough and
c           S. R. C. Malin, Report No. 71/1, Institute of Geological
c           Sciences, U.K.

c     NOSC modifications:

c        Does not do interpolation or time rate of change; calculations
c        are always for the model date, chosen to to be for the latest
c        epoch (1985).

c        Uses program INIT_GM to create a set of data statements to
c        initialize the coefficients to the most recent model: IGRF85.

c        Uses normalization coefficients generated by LEGENDRE to
c        modify the magnetic field coefficients for use by NEWMAG85.

c        The altitude range is -1 to 600 km.

c        Timings (results of 6643 calls to the subroutines):
c           NEWMAG    12 seconds; original program by McIlwain; 9 terms
c           GEOMAG85  52 seconds; using nested   do loops; 10 terms
c           GEOMAG85  39 seconds; using expanded do loops; 10 terms
c           NEWMAG85  11 seconds; 10 terms

c  Parameters passed:
c     spheroid      [i] =1: use spheroidal earth
c     alt           [r] altitude at which to calculate field; km
c     rlng          [r] West longitude; radians
c     rclt          [r] co-latitude; radians

c  Parameters returned:
c     bmf           [r] declination of the geomagnetic field; radians
c     dip           [r] dip angle; radians
c     b             [r] total field;            Webers/m**2
c     br            [r] radial component;       Webers/m**2
c     bp            [r] longitudinal component; Webers/m**2
c     bt            [r] latitudinal component;  Webers/m**2

c  Common blocks referenced:

c  Functions and subroutines referenced:
c     abs
c     acos
c     cos
c     min
c     sin
c     sqrt

c  References:

c  Change History:

c*******************!***************************************************

      integer       spheroid

      dimension     gh(0:10,0:10)

      data          nmax/10/
      data          gh/0.,
     & -2.9877000E+04,-3.1095000E+03, 3.2500000E+03, 4.0993750E+03,
     & -1.6931250E+03, 7.5075000E+02, 2.0109375E+03, 1.0557422E+03,
     &  4.7480469E+02,-7.2170313E+02, 5.4970000E+03,-1.9030000E+03,
     &  5.2740947E+03,-6.7605917E+03, 4.3165090E+03, 3.6193029E+03,
     &  1.2287031E+03,-2.1636458E+03, 4.0218750E+02, 1.2740347E+03,
     & -9.7314430E+02,-3.7949233E+03,-2.6760185E+02, 1.4644490E+03,
     &  2.4089956E+03, 1.4204622E+03, 1.9443589E+03, 7.4721161E+02,
     &  5.7921620E+01, 0.0000000E+00, 1.0865004E+02, 4.2138384E+02,
     & -9.5530100E+02, 5.4996364E+02,-2.3400855E+02, 6.6012546E+02,
     & -8.9104293E+02,-4.4238399E+02,-1.8530848E+03, 4.9148124E+02,
     & -4.5561531E+02,-9.9579408E+02,-8.2640170E+02, 1.2894187E+03,
     & -9.7827974E+02, 1.4223220E+02,-2.2037397E+02, 1.2497719E+02,
     & -3.5718332E+02, 2.1827448E+01,-7.4093585E+01,-2.4062598E+02,
     &  5.0738165E+02,-2.3374170E+02, 4.7782932E+02, 1.1374115E+03,
     & -7.2946296E+02,-1.6638974E+02, 6.6648272E+01,-3.3674916E+01,
     &  3.9555835E+01, 2.4697862E+01, 2.9661173E+01,-1.0107284E+02,
     &  3.6957808E+02,-3.0245000E+02, 1.3449809E+03, 6.8743468E+02,
     & -2.7284310E+02,-9.3072552E+00, 1.3433866E+01,-6.8512716E+01,
     &  2.1796421E+01, 2.7460910E+01,-1.7397931E+01, 1.2396026E+02,
     & -2.9085075E+03,-7.5298106E+02,-2.0478385E+01, 2.8402541E+02,
     &  1.0496591E+02,-5.0858317E+01,-3.8835591E+00, 0.0000000E+00,
     &  1.0027306E+01, 5.2734674E+01, 2.0043185E+01, 4.6921875E+02,
     & -1.1777297E+03, 2.0709787E+02,-6.6840549E+02, 1.6313645E+02,
     &  8.2382729E+01,-4.0109226E+01,-6.2670665E+00,-3.7602399E+00,
     &  5.1679555E+00, 1.6365192E+01,-2.6754728E+03, 1.7384007E+03,
     &  7.4684556E+02,-2.8187869E+02,-2.0214569E+02, 1.5658138E+02,
     &  7.5335249E+01,-1.5503866E+01, 1.2180988E+00,-3.0452470E+00,
     &  7.9643543E+00, 2.4328607E+02, 0.0000000E+00, 4.9584102E+02,
     &  7.0122510E+02,-2.9566246E+02, 0.0000000E+00,-2.0043185E+01,
     &  3.2730385E+01, 0.0000000E+00,-3.5617675E+00, 0.0000000E+00/

c     A2 and B2 are the squares of the semi-major and semi-minor axes
c     of the reference spheroid used to transform to geodetic from
c     geocentric coordinates.

      data          a2 /40680925./, b2 /40408588./

      REDUCE(arg)=SIGN(MIN(ABS(arg),1.),arg)


c     Terms in THETA (co-latitude)
      st=ABS(SIN(rclt))
      if (st .lt. 1.e-6) then
          st=1.e-6
          nmax=5
      end if
      ct=SQRT(1.-st*st)
      if (rclt .gt. 1.570796327e0) ct=-ct

c     Terms in PHI (East longitude)
      phi=rlng
      if (phi .lt. 0.) then
         phi=ABS(phi)
      else
     &if (phi .gt. 0.) then
         phi=6.2831853e0-phi
      end if
      sp 1=SIN(phi)
      cp 1=COS(phi)

      if (spheroid .eq. 0) then

c        No spheroidal correction
         r=alt+6371.2
         cd=1.
         sd=0.
      else

c        Spheroidal correction
         aa=a2*st*st
         bb=b2*ct*ct
         cc=SQRT(aa+bb)
         r=SQRT(alt*(alt+2.*cc)+(a2*aa+b2*bb)/(aa+bb))
         cd=(alt+cc)/r
         sd=(a2-b2)/cc*st*ct/r
         aa=ct
         ct=ct*cd-st*sd
         st=st*cd+aa*sd
      end if
      ar=6371.2/r

c     n= 1: First order terms
      aor=ar*ar*ar

      c 1=gh( 1, 1)*cp 1+gh( 0, 1)*sp 1

      Q 1 0=ct
      Q 1 1=st

      dQ 1 0=-st
      dQ 1 1= ct

      br= aor*( Q 1 0*gh( 1, 0)+ Q 1 1*c 1)* 2.

      bt=-aor*(dQ 1 0*gh( 1, 0)+dQ 1 1*c 1)

      bp= aor*  Q 1 1*(gh( 1, 1)*sp 1-gh( 0, 1)*cp 1)

      if (nmax .lt. 2) go to 260

      n= 2
      aor=aor*ar

      sp 2=(sp 1+sp 1)*cp 1
      cp 2=(cp 1+sp 1)*(cp 1-sp 1)

      c 1=gh( 2, 1)*cp 1+gh( 0, 2)*sp 1
      c 2=gh( 2, 2)*cp 2+gh( 1, 2)*sp 2

      Q 2 0=ct*ct-.33333333
      Q 2 1=st*ct
      Q 2 2=st*st

      dQ 2 0=-st*ct-st*ct
      dQ 2 1= ct*ct-st*st
      dQ 2 2=-dQ 2 0

      br=br+aor*( Q 2 0*gh( 2, 0)+ Q 2 1*c 1+ Q 2 2*c 2)* 3.

      bt=bt-aor*(dQ 2 0*gh( 2, 0)+dQ 2 1*c 1+dQ 2 2*c 2)

      bp=bp+aor*( Q 2 1*(gh( 2, 1)*sp 1-gh( 0, 2)*cp 1)
     &          + Q 2 2*(gh( 2, 2)*sp 2-gh( 1, 2)*cp 2)* 2.)

      if (nmax .lt.  3) go to 260

      n= 3
      aor=aor*ar

      sp 3=sp 1*cp 2+cp 1*sp 2
      cp 3=cp 1*cp 2-sp 1*sp 2

      c 1=gh( 3, 1)*cp 1+gh( 0, 3)*sp 1
      c 2=gh( 3, 2)*cp 2+gh( 1, 3)*sp 2
      c 3=gh( 3, 3)*cp 3+gh( 2, 3)*sp 3

      Q 3 0=ct*Q 2 0-.26666667*Q 1 0
      Q 3 1=ct*Q 2 1-.20000000*Q 1 1
      Q 3 2=ct*Q 2 2
      Q 3 3=st*Q 2 2

      dQ 3 0=ct*dQ 2 0-st*Q 2 0-.26666667*dQ 1 0
      dQ 3 1=ct*dQ 2 1-st*Q 2 1-.20000000*dQ 1 1
      dQ 3 2=ct*dQ 2 2-st*Q 2 2
      dQ 3 3= 3.*Q 3 2

      br=br+aor*( Q 3 0*gh( 3, 0)+ Q 3 1*c 1+ Q 3 2*c 2+ Q 3 3*c 3)* 4.

      bt=bt-aor*(dQ 3 0*gh( 3, 0)+dQ 3 1*c 1+dQ 3 2*c 2+dQ 3 3*c 3)

      bp=bp+aor*( Q 3 1*(gh( 3, 1)*sp 1-gh( 0, 3)*cp 1)
     &          + Q 3 2*(gh( 3, 2)*sp 2-gh( 1, 3)*cp 2)* 2.
     &          + Q 3 3*(gh( 3, 3)*sp 3-gh( 2, 3)*cp 3)* 3.)

      if (nmax .lt.  4) go to 260

      n= 4
      aor=aor*ar

      sp 4=sp 1*cp 3+cp 1*sp 3
      cp 4=cp 1*cp 3-sp 1*sp 3

      c 1=gh( 4, 1)*cp 1+gh( 0, 4)*sp 1
      c 2=gh( 4, 2)*cp 2+gh( 1, 4)*sp 2
      c 3=gh( 4, 3)*cp 3+gh( 2, 4)*sp 3
      c 4=gh( 4, 4)*cp 4+gh( 3, 4)*sp 4

      Q 4 0=ct*Q 3 0-.25714286*Q 2 0
      Q 4 1=ct*Q 3 1-.22857143*Q 2 1
      Q 4 2=ct*Q 3 2-.14285714*Q 2 2
      Q 4 3=ct*Q 3 3
      Q 4 4=st*Q 3 3

      dQ 4 0=ct*dQ 3 0-st*Q 3 0-.25714286*dQ 2 0
      dQ 4 1=ct*dQ 3 1-st*Q 3 1-.22857143*dQ 2 1
      dQ 4 2=ct*dQ 3 2-st*Q 3 2-.14285714*dQ 2 2
      dQ 4 3=ct*dQ 3 3-st*Q 3 3
      dQ 4 4= 4.*Q 4 3

      br=br+aor*( Q 4 0*gh( 4, 0)+ Q 4 1*c 1+ Q 4 2*c 2+ Q 4 3*c 3
     &                           + Q 4 4*c 4)* 5.

      bt=bt-aor*(dQ 4 0*gh( 4, 0)+dQ 4 1*c 1+dQ 4 2*c 2+dQ 4 3*c 3
     &                           +dQ 4 4*c 4)

      bp=bp+aor*( Q 4 1*(gh( 4, 1)*sp 1-gh( 0, 4)*cp 1)
     &          + Q 4 2*(gh( 4, 2)*sp 2-gh( 1, 4)*cp 2)* 2.
     &          + Q 4 3*(gh( 4, 3)*sp 3-gh( 2, 4)*cp 3)* 3.
     &          + Q 4 4*(gh( 4, 4)*sp 4-gh( 3, 4)*cp 4)* 4.)

      if (nmax .lt.  5) go to 260

      n= 5
      aor=aor*ar

      sp 5=sp 1*cp 4+cp 1*sp 4
      cp 5=cp 1*cp 4-sp 1*sp 4

      c 1=gh( 5, 1)*cp 1+gh( 0, 5)*sp 1
      c 2=gh( 5, 2)*cp 2+gh( 1, 5)*sp 2
      c 3=gh( 5, 3)*cp 3+gh( 2, 5)*sp 3
      c 4=gh( 5, 4)*cp 4+gh( 3, 5)*sp 4
      c 5=gh( 5, 5)*cp 5+gh( 4, 5)*sp 5

      Q 5 0=ct*Q 4 0-.25396825*Q 3 0
      Q 5 1=ct*Q 4 1-.23809524*Q 3 1
      Q 5 2=ct*Q 4 2-.19047619*Q 3 2
      Q 5 3=ct*Q 4 3-.11111111*Q 3 3
      Q 5 4=ct*Q 4 4
      Q 5 5=st*Q 4 4

      dQ 5 0=ct*dQ 4 0-st*Q 4 0-.25396825*dQ 3 0
      dQ 5 1=ct*dQ 4 1-st*Q 4 1-.23809524*dQ 3 1
      dQ 5 2=ct*dQ 4 2-st*Q 4 2-.19047619*dQ 3 2
      dQ 5 3=ct*dQ 4 3-st*Q 4 3-.11111111*dQ 3 3
      dQ 5 4=ct*dQ 4 4-st*Q 4 4
      dQ 5 5= 5.*Q 5 4

      br=br+aor*( Q 5 0*gh( 5, 0)+ Q 5 1*c 1+ Q 5 2*c 2+ Q 5 3*c 3
     &                           + Q 5 4*c 4+ Q 5 5*c 5)* 6.

      bt=bt-aor*(dQ 5 0*gh( 5, 0)+dQ 5 1*c 1+dQ 5 2*c 2+dQ 5 3*c 3
     &                           +dQ 5 4*c 4+dQ 5 5*c 5)

      bp=bp+aor*( Q 5 1*(gh( 5, 1)*sp 1-gh( 0, 5)*cp 1)
     &          + Q 5 2*(gh( 5, 2)*sp 2-gh( 1, 5)*cp 2)* 2.
     &          + Q 5 3*(gh( 5, 3)*sp 3-gh( 2, 5)*cp 3)* 3.
     &          + Q 5 4*(gh( 5, 4)*sp 4-gh( 3, 5)*cp 4)* 4.
     &          + Q 5 5*(gh( 5, 5)*sp 5-gh( 4, 5)*cp 5)* 5.)

      if (nmax .lt.  6) go to 260

      n= 6
      aor=aor*ar

      sp 6=sp 1*cp 5+cp 1*sp 5
      cp 6=cp 1*cp 5-sp 1*sp 5

      c 1=gh( 6, 1)*cp 1+gh( 0, 6)*sp 1
      c 2=gh( 6, 2)*cp 2+gh( 1, 6)*sp 2
      c 3=gh( 6, 3)*cp 3+gh( 2, 6)*sp 3
      c 4=gh( 6, 4)*cp 4+gh( 3, 6)*sp 4
      c 5=gh( 6, 5)*cp 5+gh( 4, 6)*sp 5
      c 6=gh( 6, 6)*cp 6+gh( 5, 6)*sp 6

      Q 6 0=ct*Q 5 0-.25252525*Q 4 0
      Q 6 1=ct*Q 5 1-.24242424*Q 4 1
      Q 6 2=ct*Q 5 2-.21212121*Q 4 2
      Q 6 3=ct*Q 5 3-.16161616*Q 4 3
      Q 6 4=ct*Q 5 4-.09090909*Q 4 4
      Q 6 5=ct*Q 5 5
      Q 6 6=st*Q 5 5

      dQ 6 0=ct*dQ 5 0-st*Q 5 0-.25252525*dQ 4 0
      dQ 6 1=ct*dQ 5 1-st*Q 5 1-.24242424*dQ 4 1
      dQ 6 2=ct*dQ 5 2-st*Q 5 2-.21212121*dQ 4 2
      dQ 6 3=ct*dQ 5 3-st*Q 5 3-.16161616*dQ 4 3
      dQ 6 4=ct*dQ 5 4-st*Q 5 4-.09090909*dQ 4 4
      dQ 6 5=ct*dQ 5 5-st*Q 5 5
      dQ 6 6= 6.*Q 6 5

      br=br+aor*( Q 6 0*gh( 6, 0)+ Q 6 1*c 1+ Q 6 2*c 2+ Q 6 3*c 3
     &                           + Q 6 4*c 4+ Q 6 5*c 5+ Q 6 6*c 6)* 7.

      bt=bt-aor*(dQ 6 0*gh( 6, 0)+dQ 6 1*c 1+dQ 6 2*c 2+dQ 6 3*c 3
     &                           +dQ 6 4*c 4+dQ 6 5*c 5+dQ 6 6*c 6)

      bp=bp+aor*( Q 6 1*(gh( 6, 1)*sp 1-gh( 0, 6)*cp 1)
     &          + Q 6 2*(gh( 6, 2)*sp 2-gh( 1, 6)*cp 2)* 2.
     &          + Q 6 3*(gh( 6, 3)*sp 3-gh( 2, 6)*cp 3)* 3.
     &          + Q 6 4*(gh( 6, 4)*sp 4-gh( 3, 6)*cp 4)* 4.
     &          + Q 6 5*(gh( 6, 5)*sp 5-gh( 4, 6)*cp 5)* 5.
     &          + Q 6 6*(gh( 6, 6)*sp 6-gh( 5, 6)*cp 6)* 6.)

      if (nmax .lt.  7) go to 260

      n= 7
      aor=aor*ar

      sp 7=sp 1*cp 6+cp 1*sp 6
      cp 7=cp 1*cp 6-sp 1*sp 6

      c 1=gh( 7, 1)*cp 1+gh( 0, 7)*sp 1
      c 2=gh( 7, 2)*cp 2+gh( 1, 7)*sp 2
      c 3=gh( 7, 3)*cp 3+gh( 2, 7)*sp 3
      c 4=gh( 7, 4)*cp 4+gh( 3, 7)*sp 4
      c 5=gh( 7, 5)*cp 5+gh( 4, 7)*sp 5
      c 6=gh( 7, 6)*cp 6+gh( 5, 7)*sp 6
      c 7=gh( 7, 7)*cp 7+gh( 6, 7)*sp 7

      Q 7 0=ct*Q 6 0-.25174825*Q 5 0
      Q 7 1=ct*Q 6 1-.24475524*Q 5 1
      Q 7 2=ct*Q 6 2-.22377622*Q 5 2
      Q 7 3=ct*Q 6 3-.18881119*Q 5 3
      Q 7 4=ct*Q 6 4-.13986014*Q 5 4
      Q 7 5=ct*Q 6 5-.07692308*Q 5 5
      Q 7 6=ct*Q 6 6
      Q 7 7=st*Q 6 6

      dQ 7 0=ct*dQ 6 0-st*Q 6 0-.25174825*dQ 5 0
      dQ 7 1=ct*dQ 6 1-st*Q 6 1-.24475524*dQ 5 1
      dQ 7 2=ct*dQ 6 2-st*Q 6 2-.22377622*dQ 5 2
      dQ 7 3=ct*dQ 6 3-st*Q 6 3-.18881119*dQ 5 3
      dQ 7 4=ct*dQ 6 4-st*Q 6 4-.13986014*dQ 5 4
      dQ 7 5=ct*dQ 6 5-st*Q 6 5-.07692308*dQ 5 5
      dQ 7 6=ct*dQ 6 6-st*Q 6 6
      dQ 7 7= 7.*Q 7 6

      br=br+aor*( Q 7 0*gh( 7, 0)+ Q 7 1*c 1+ Q 7 2*c 2+ Q 7 3*c 3
     &                           + Q 7 4*c 4+ Q 7 5*c 5+ Q 7 6*c 6
     &                           + Q 7 7*c 7)* 8.

      bt=bt-aor*(dQ 7 0*gh( 7, 0)+dQ 7 1*c 1+dQ 7 2*c 2+dQ 7 3*c 3
     &                           +dQ 7 4*c 4+dQ 7 5*c 5+dQ 7 6*c 6
     &                           +dQ 7 7*c 7)

      bp=bp+aor*( Q 7 1*(gh( 7, 1)*sp 1-gh( 0, 7)*cp 1)
     &          + Q 7 2*(gh( 7, 2)*sp 2-gh( 1, 7)*cp 2)* 2.
     &          + Q 7 3*(gh( 7, 3)*sp 3-gh( 2, 7)*cp 3)* 3.
     &          + Q 7 4*(gh( 7, 4)*sp 4-gh( 3, 7)*cp 4)* 4.
     &          + Q 7 5*(gh( 7, 5)*sp 5-gh( 4, 7)*cp 5)* 5.
     &          + Q 7 6*(gh( 7, 6)*sp 6-gh( 5, 7)*cp 6)* 6.
     &          + Q 7 7*(gh( 7, 7)*sp 7-gh( 6, 7)*cp 7)* 7.)

      if (nmax .lt.  8) go to 260

      n= 8
      aor=aor*ar

      sp 8=sp 1*cp 7+cp 1*sp 7
      cp 8=cp 1*cp 7-sp 1*sp 7

      c 1=gh( 8, 1)*cp 1+gh( 0, 8)*sp 1
      c 2=gh( 8, 2)*cp 2+gh( 1, 8)*sp 2
      c 3=gh( 8, 3)*cp 3+gh( 2, 8)*sp 3
      c 4=gh( 8, 4)*cp 4+gh( 3, 8)*sp 4
      c 5=gh( 8, 5)*cp 5+gh( 4, 8)*sp 5
      c 6=gh( 8, 6)*cp 6+gh( 5, 8)*sp 6
      c 7=gh( 8, 7)*cp 7+gh( 6, 8)*sp 7
      c 8=gh( 8, 8)*cp 8+gh( 7, 8)*sp 8

      Q 8 0=ct*Q 7 0-.25128205*Q 6 0
      Q 8 1=ct*Q 7 1-.24615385*Q 6 1
      Q 8 2=ct*Q 7 2-.23076923*Q 6 2
      Q 8 3=ct*Q 7 3-.20512821*Q 6 3
      Q 8 4=ct*Q 7 4-.16923077*Q 6 4
      Q 8 5=ct*Q 7 5-.12307692*Q 6 5
      Q 8 6=ct*Q 7 6-.06666667*Q 6 6
      Q 8 7=ct*Q 7 7
      Q 8 8=st*Q 7 7

      dQ 8 0=ct*dQ 7 0-st*Q 7 0-.25128205*dQ 6 0
      dQ 8 1=ct*dQ 7 1-st*Q 7 1-.24615385*dQ 6 1
      dQ 8 2=ct*dQ 7 2-st*Q 7 2-.23076923*dQ 6 2
      dQ 8 3=ct*dQ 7 3-st*Q 7 3-.20512821*dQ 6 3
      dQ 8 4=ct*dQ 7 4-st*Q 7 4-.16923077*dQ 6 4
      dQ 8 5=ct*dQ 7 5-st*Q 7 5-.12307692*dQ 6 5
      dQ 8 6=ct*dQ 7 6-st*Q 7 6-.06666667*dQ 6 6
      dQ 8 7=ct*dQ 7 7-st*Q 7 7
      dQ 8 8= 8.*Q 8 7

      br=br+aor*( Q 8 0*gh( 8, 0)+ Q 8 1*c 1+ Q 8 2*c 2+ Q 8 3*c 3
     &                           + Q 8 4*c 4+ Q 8 5*c 5+ Q 8 6*c 6
     &                           + Q 8 7*c 7+ Q 8 8*c 8)* 9.

      bt=bt-aor*(dQ 8 0*gh( 8, 0)+dQ 8 1*c 1+dQ 8 2*c 2+dQ 8 3*c 3
     &                           +dQ 8 4*c 4+dQ 8 5*c 5+dQ 8 6*c 6
     &                           +dQ 8 7*c 7+dQ 8 8*c 8)

      bp=bp+aor*( Q 8 1*(gh( 8, 1)*sp 1-gh( 0, 8)*cp 1)
     &          + Q 8 2*(gh( 8, 2)*sp 2-gh( 1, 8)*cp 2)* 2.
     &          + Q 8 3*(gh( 8, 3)*sp 3-gh( 2, 8)*cp 3)* 3.
     &          + Q 8 4*(gh( 8, 4)*sp 4-gh( 3, 8)*cp 4)* 4.
     &          + Q 8 5*(gh( 8, 5)*sp 5-gh( 4, 8)*cp 5)* 5.
     &          + Q 8 6*(gh( 8, 6)*sp 6-gh( 5, 8)*cp 6)* 6.
     &          + Q 8 7*(gh( 8, 7)*sp 7-gh( 6, 8)*cp 7)* 7.
     &          + Q 8 8*(gh( 8, 8)*sp 8-gh( 7, 8)*cp 8)* 8.)

      if (nmax .lt.  9) go to 260

      n= 9
      aor=aor*ar

      sp 9=sp 1*cp 8+cp 1*sp 8
      cp 9=cp 1*cp 8-sp 1*sp 8

      c 1=gh( 9, 1)*cp 1+gh( 0, 9)*sp 1
      c 2=gh( 9, 2)*cp 2+gh( 1, 9)*sp 2
      c 3=gh( 9, 3)*cp 3+gh( 2, 9)*sp 3
      c 4=gh( 9, 4)*cp 4+gh( 3, 9)*sp 4
      c 5=gh( 9, 5)*cp 5+gh( 4, 9)*sp 5
      c 6=gh( 9, 6)*cp 6+gh( 5, 9)*sp 6
      c 7=gh( 9, 7)*cp 7+gh( 6, 9)*sp 7
      c 8=gh( 9, 8)*cp 8+gh( 7, 9)*sp 8
      c 9=gh( 9, 9)*cp 9+gh( 8, 9)*sp 9

      Q 9 0=ct*Q 8 0-.25098039*Q 7 0
      Q 9 1=ct*Q 8 1-.24705882*Q 7 1
      Q 9 2=ct*Q 8 2-.23529412*Q 7 2
      Q 9 3=ct*Q 8 3-.21568627*Q 7 3
      Q 9 4=ct*Q 8 4-.18823529*Q 7 4
      Q 9 5=ct*Q 8 5-.15294118*Q 7 5
      Q 9 6=ct*Q 8 6-.10980392*Q 7 6
      Q 9 7=ct*Q 8 7-.05882353*Q 7 7
      Q 9 8=ct*Q 8 8
      Q 9 9=st*Q 8 8

      dQ 9 0=ct*dQ 8 0-st*Q 8 0-.25098039*dQ 7 0
      dQ 9 1=ct*dQ 8 1-st*Q 8 1-.24705882*dQ 7 1
      dQ 9 2=ct*dQ 8 2-st*Q 8 2-.23529412*dQ 7 2
      dQ 9 3=ct*dQ 8 3-st*Q 8 3-.21568627*dQ 7 3
      dQ 9 4=ct*dQ 8 4-st*Q 8 4-.18823529*dQ 7 4
      dQ 9 5=ct*dQ 8 5-st*Q 8 5-.15294118*dQ 7 5
      dQ 9 6=ct*dQ 8 6-st*Q 8 6-.10980392*dQ 7 6
      dQ 9 7=ct*dQ 8 7-st*Q 8 7-.05882353*dQ 7 7
      dQ 9 8=ct*dQ 8 8-st*Q 8 8
      dQ 9 9= 9.*Q 9 8

      br=br+aor*( Q 9 0*gh( 9, 0)+ Q 9 1*c 1+ Q 9 2*c 2+ Q 9 3*c 3
     &                           + Q 9 4*c 4+ Q 9 5*c 5+ Q 9 6*c 6
     &                           + Q 9 7*c 7+ Q 9 8*c 8+ Q 9 9*c 9)*10.

      bt=bt-aor*(dQ 9 0*gh( 9, 0)+dQ 9 1*c 1+dQ 9 2*c 2+dQ 9 3*c 3
     &                           +dQ 9 4*c 4+dQ 9 5*c 5+dQ 9 6*c 6
     &                           +dQ 9 7*c 7+dQ 9 8*c 8+dQ 9 9*c 9)

      bp=bp+aor*( Q 9 1*(gh( 9, 1)*sp 1-gh( 0, 9)*cp 1)
     &          + Q 9 2*(gh( 9, 2)*sp 2-gh( 1, 9)*cp 2)* 2.
     &          + Q 9 3*(gh( 9, 3)*sp 3-gh( 2, 9)*cp 3)* 3.
     &          + Q 9 4*(gh( 9, 4)*sp 4-gh( 3, 9)*cp 4)* 4.
     &          + Q 9 5*(gh( 9, 5)*sp 5-gh( 4, 9)*cp 5)* 5.
     &          + Q 9 6*(gh( 9, 6)*sp 6-gh( 5, 9)*cp 6)* 6.
     &          + Q 9 7*(gh( 9, 7)*sp 7-gh( 6, 9)*cp 7)* 7.
     &          + Q 9 8*(gh( 9, 8)*sp 8-gh( 7, 9)*cp 8)* 8.
     &          + Q 9 9*(gh( 9, 9)*sp 9-gh( 8, 9)*cp 9)* 9.)

      if (nmax .lt. 10) go to 260

      n=10
      aor=aor*ar

      sp10=sp 1*cp 9+cp 1*sp 9
      cp10=cp 1*cp 9-sp 1*sp 9

      c 1=gh(10, 1)*cp 1+gh( 0,10)*sp 1
      c 2=gh(10, 2)*cp 2+gh( 1,10)*sp 2
      c 3=gh(10, 3)*cp 3+gh( 2,10)*sp 3
      c 4=gh(10, 4)*cp 4+gh( 3,10)*sp 4
      c 5=gh(10, 5)*cp 5+gh( 4,10)*sp 5
      c 6=gh(10, 6)*cp 6+gh( 5,10)*sp 6
      c 7=gh(10, 7)*cp 7+gh( 6,10)*sp 7
      c 8=gh(10, 8)*cp 8+gh( 7,10)*sp 8
      c 9=gh(10, 9)*cp 9+gh( 8,10)*sp 9
      c10=gh(10,10)*cp10+gh( 9,10)*sp10

      Q10 0=ct*Q 9 0-.25077399*Q 8 0
      Q10 1=ct*Q 9 1-.24767802*Q 8 1
      Q10 2=ct*Q 9 2-.23839009*Q 8 2
      Q10 3=ct*Q 9 3-.22291022*Q 8 3
      Q10 4=ct*Q 9 4-.20123839*Q 8 4
      Q10 5=ct*Q 9 5-.17337461*Q 8 5
      Q10 6=ct*Q 9 6-.13931889*Q 8 6
      Q10 7=ct*Q 9 7-.09907121*Q 8 7
      Q10 8=ct*Q 9 8-.05263158*Q 8 8
      Q10 9=ct*Q 9 9
      Q1010=st*Q 9 9

      dQ10 0=ct*dQ 9 0-st*Q 9 0-.25077399*dQ 8 0
      dQ10 1=ct*dQ 9 1-st*Q 9 1-.24767802*dQ 8 1
      dQ10 2=ct*dQ 9 2-st*Q 9 2-.23839009*dQ 8 2
      dQ10 3=ct*dQ 9 3-st*Q 9 3-.22291022*dQ 8 3
      dQ10 4=ct*dQ 9 4-st*Q 9 4-.20123839*dQ 8 4
      dQ10 5=ct*dQ 9 5-st*Q 9 5-.17337461*dQ 8 5
      dQ10 6=ct*dQ 9 6-st*Q 9 6-.13931889*dQ 8 6
      dQ10 7=ct*dQ 9 7-st*Q 9 7-.09907121*dQ 8 7
      dQ10 8=ct*dQ 9 8-st*Q 9 8-.05263158*dQ 8 8
      dQ10 9=ct*dQ 9 9-st*Q 9 9
      dQ1010=10.*Q10 9

      br=br+aor*( Q10 0*gh(10, 0)+ Q10 1*c 1+ Q10 2*c 2+ Q10 3*c 3
     &                           + Q10 4*c 4+ Q10 5*c 5+ Q10 6*c 6
     &                           + Q10 7*c 7+ Q10 8*c 8+ Q10 9*c 9
     &                           + Q1010*c10)*11.

      bt=bt-aor*(dQ10 0*gh(10, 0)+dQ10 1*c 1+dQ10 2*c 2+dQ10 3*c 3
     &                           +dQ10 4*c 4+dQ10 5*c 5+dQ10 6*c 6
     &                           +dQ10 7*c 7+dQ10 8*c 8+dQ10 9*c 9
     &                           +dQ1010*c10)

      bp=bp+aor*( Q10 1*(gh(10, 1)*sp 1-gh( 0,10)*cp 1)
     &          + Q10 2*(gh(10, 2)*sp 2-gh( 1,10)*cp 2)* 2.
     &          + Q10 3*(gh(10, 3)*sp 3-gh( 2,10)*cp 3)* 3.
     &          + Q10 4*(gh(10, 4)*sp 4-gh( 3,10)*cp 4)* 4.
     &          + Q10 5*(gh(10, 5)*sp 5-gh( 4,10)*cp 5)* 5.
     &          + Q10 6*(gh(10, 6)*sp 6-gh( 5,10)*cp 6)* 6.
     &          + Q10 7*(gh(10, 7)*sp 7-gh( 6,10)*cp 7)* 7.
     &          + Q10 8*(gh(10, 8)*sp 8-gh( 7,10)*cp 8)* 8.
     &          + Q10 9*(gh(10, 9)*sp 9-gh( 8,10)*cp 9)* 9.
     &          + Q1010*(gh(10,10)*sp10-gh( 9,10)*cp10)*10.)

260   bp=bp*1.e-5/st
      bt=bt*1.e-5
      br=br*1.e-5

      tmp=bt
      bt=cd*bt+sd*br
      br=cd*br-sd*tmp

      b=SQRT(br*br+bt*bt+bp*bp)
      bh=SQRT(bt*bt+bp*bp)
      bmf=3.141592654e0-ACOS(REDUCE(bt/bh))
      if (bp .lt. 0.) bmf=-bmf
      dip=ACOS(REDUCE(bh/b))
      if (br .gt. 0.) dip=-dip

      RETURN
      END      !  NEWMAG
